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Abstract 

We show that discrete schemes developed for lattice hydrodynamics pro- 
vide an elegant and physically transparent way of deriving Laplacians with 
isotropic discretisation error. Isotropy is guaranteed whenever the Laplacian 
weights follow from the discrete Maxwell-Boltzmann equilibrium since these 
are, by construction, isotropic on the lattice. We also point out that stencils 
using as few as 15 points in three dimensions, generate isotropic Laplacians. 
These computationally efficient Laplacians can be used in cell-dynamical and 
hybrid lattice Boltzmann simulations, in favor of popular anisotropic Lapla- 
cians, which make use of larger stencils. The method can be extended to 
provide discretisations of higher order and for other differential operators, 
such the gradient, divergence and curl. 



1. Introduction 

Isotropy, an essential property of the Laplacian operator, is a desirable 
feature in any of its discrete representations. However, most commonly used 
discretisations of the Laplacian suffer from error terms that are anisotropic. 
These infect numerical solutions with anisotropies when, in fact, physical 
solutions are required to be isotropic. Two recent studies have addressed this 
issue, by providing algebraic methods for constructing isotropic Laplacians 
[H, lij. In this paper, we show that lattices and associated weights used in 
lattice hydrodynamic simulations, naturally provide discrete Laplacians with 
isotropic discretisation error. These lattices and weights have been known for 
a long time in the lattice Boltzmann literature, jsl, |J, but their connection to 
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Figure 1: Ordering of points on a cubic unit cell according to "energy shells" , 
as explained in the text. Here 69,61,62,63 represent the energy shells corre- 
sponding to 6 = 0, 1, 2, 3, marked by o, •, x and ■ respectively. For clarity, 
only one pair is highlighted. 



isotropic Laplacians has not been recognized. Here, we make this connection 
explicit and also show how discretisations of other differential operators like 
the gradient and divergence follow from it. 

Below, we first provide our main results for the isotropically discretised 
Laplacian. We then provide a derivation of our result and make explicit its 
connection with the discrete Maxwell-Boltzmann equilibrium on the lattice. 
We conclude with a comparison of the isotropy properties of our Laplacians 
with those proposed earlier. 

Consider the cubic cell in Fig. [1] with center 60, the 6 vectors c^^^ that 
point to the face centers 61, the 12 vectors that point to the edge centers 
62 and the 8 vectors cf that point to the vertices 63. The face center vectors 
point along the Cartesian axes and so have the form (±1,0,0) with two 
vanishing Cartesian components, the edge center vectors are confined to the 
Cartesian planes and so have form (±1, ±1, 0) with one vanishing Cartesian 
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component, while the vertex vectors have the form (±1, ±1, ±1). (With this 
choice, the discretisation step is set to unity.) This suggests an "energy shell" 
classification, where the energy is identified with the squared modulus of the 
vector. Then, the 27 points of the cubic cell lie on the "energy shells" with 
energies ej = 0, 1, 2, 3. There is 1 null vector c° on the energy shell e = 0, 6 
vectors cj on the energy shell e = 1, 12 vectors cf on the energy shell e = 2 
and 8 vectors on the energy shell e = 3. In the cell dynamics literature 
these are known as the nearest neighbors - NN, next nearest neighbors - 
NNN, and next next nearest neighbors - NNNN, respectively. 

For a function ip{r) defined on such a cubic grid, we use the convenient 
shorthand iplr + c^) = tp^''^ for its value on the z-th point of the j-th energy 
shell. Then, our main results for the isotropic Laplacian L(r) = V'^tp{r) are 
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The first of these is a two-dimensional Laplacian which uses the 4 face centers 
and 4 edge centers on a Cartesian plane, while the remaining are all three 
dimensional Laplacians. The subscript on each of these Laplacians indicates 
that they derive from a DnQm lattice hydrodynamic model, as we explain 
in detail below. 



2. Laplacians from lattice equilibria 

The lattice formulation of kinetic theory provides a computationally ef- 
ficient algorithm for solving the Navier-Stokes and related equations A 
central quantity in lattice kinetic theory is the discrete form of the Maxwell- 
Boltzmann velocity distribution. There are several routes by which these 
equilibria can be obtained jsl, 0]. Here, we focus on the lattice generalization 
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of the earliest derivation by Maxwell [5| , which appeals only to factorisabil- 
ity and isotropy. Briefly, Maxwell argued that the distribution function /(c) 
must be isotropic in velocity space and so must be function of = c ■ c and 
it must be factorisable so that /(c) = f{cx)f{cy)f{cz). The only function 
that satisfles both these requirements is the Gaussian /(c) ~ exp(— c^/2/cbT), 
whose variance fc^T is fixed by requiring consistency with thermodynamics, 
and whose prefactor is fixed by normalization. In the above, the mass has 
been set to unity. Recent works 0, S] show that a cubic lattice containing 
27 velocities, made out of the direct product of the velocity set {—1,0,1}, 
has both the isotropy and factorisability features required by the continuum 
Maxwellian. This provides a succinct derivation of the lattice Maxwellian. 
The 27 velocity lattice and the weights so generated are called D3Q27 in the 
lattice hydrodynamics literature. Projections of this lattice to smaller veloc- 
ity set are possible. In such projections, factorisability is lost, but isotropy is 
still maintained. It is well-known that isotropy of the lattice Maxwellians is 
a necessary condition for Navier-Stokes hydrodjTiamics on the lattice. It is 
this feature that we exploit to generate the isotropic Laplacians listed above. 

We begin with a DriQm lattice hydrodynamic model in n dimensions with 
m velocities. These have discrete weights wl which correspond to the energy 
shells Cj. For a scalar function ip{r) defined on such a lattice, consider the 
transform 

3 N, 
i=0 i=l 

where i labels the i-th discrete velocity in the j-th energy shell with Nj 
velocities and wl is the corresponding weight factor. As is well known, the 
necessary conditions for obtaining isotropic Navier-Stokes hydrodynamics on 
the lattice are ^ 

E^' = l (6) 

Y.''lciA, = T5^P (7) 

where Greek indices label Cartesian directions and i^^^^x ~ ^a/3^7A + ^aA^7^ + 
Sa^Spx. In the above T is a lattice-dependent constant which identifies with 
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the temperature of the moving particles. Using particles in the cubic cell 
with velocities Ci^a = {—I? 0,1}, the identity = cf^ implies T = 1/3 as 
the only value ensuring isotropy of the lattice hydrodynamics in the cubic 
cell. All weighted polynomials odd in the velocities vanish identically. In 
particular, linear, cubic and quintic polynomials are zero. On these lattices, 
the sextic polynomials are the first non-zero polynomials to break isotropy. 

Taylor expanding iplr + c^) in Eq. Eland applying the above symmetries 
of Eq. [6] - m we obtain 



^(r) + |v2v^(r) 



(9) 



This equation can be solved for Laplacian L(r) = V'^iP{t) to obtain 
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This automatically secures isotropy of the Laplacian up to leading order error, 
with an error coefficient of order 0{T). This error cannot be made zero with 
the given cubic stencil. The above expression is remarkable, because any 
lattice with suitable weights which satisfies the conditions in Eq. |6] - [HI will 
provide an expression for the discrete Laplacian operator and will ensure 
isotropy. Writing down the terms of Eq. [TO] explicitly, with ip^''^ for 'ijj(r + ci), 
we have 



L(r) 



2 
T 



12 



(3) 



[W 



1)^ 



(0) 



i=l 



1=1 



i=l 



In other words, by redefining 
have 
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1 and = for j = 1,2,3 we 
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to replace Eq. [6] which, along with Eq. [7HH1 form a set of weights, {w}, 
required to construct isotropic Laplacian operators in discrete space. These 
weights are lattice analogues of Hermite weights related to the Maxwell- 
Boltzmann equilibrium. Hence, this method of deriving isotropic Laplacians 
is an elegant and physically transparent way of calculating the operators 
compared to the methods present in the literature [ll, 0]- 
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Table 1: Energy shells and the corresponding weight factors for various 
DnQm lattice hydrodynamics models. Values of Nj for the two dimensional 
model D2Q9 are given in brackets. 

This general expression, when applied to the D2Q9, D3Q15, D3Q19, 
D3Q27 models, whose weights are listed in table [H gives the isotropic Lapla- 
cians we listed previously. The two-dimensional D2Q9 Laplacian was ob- 
tained earlier by a different argument in [l|. Interestingly, isotropic three- 
dimensional Laplacians can be achieved with just 15 or 19 velocities, i.e in 
D3Q15 and D3Q19 using Eq. [TOl and it is not necessary to use all 27 veloci- 
ties. This result was obtained earlier 10|,lll|], with weights identical to D3Q15 
and D3Q19, without realizing the connection with lattice hydrodynamics. 

3. Discussion 

Having derived the Laplacians, we next proceed to compare their isotropy 
properties with other commonly used Laplacians, listed below. 
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The suffixes CD, PK, SO, KU ,EW stand for central difference, Patra- 
Kartunnen, Sliinozaki-Oono, Kumar and 'equally weighted' respectively. Eq. 
T3] is the standard central finite-difference expression. Eq. [H] has been 
systematically derived by imposing conditions of rotational invariance and 
isotropy of the operator I, EL Eq. m is popular in the cell-dynamics 
and phase separation studies [l3|. Eq. [12] has been introduced as part 
of isotropic finite differences which describes discrete derivative operations 
without directional bias [l[ . Eq. [17] is a simple expression used in lat- 
tice Boltzmann simulations |15| which gives equal weightage to all en- 
ergy shells. In the small wave number limit, the discrete Fourier transform 
L(k) = exp(— ?k ■ r)L(r) / exp(— ik ■ r)'?/;(r) of the Laplacian operator 
corresponding to the above expressions (Eq. [2] - |H [12] - [T71) may be written 
as 
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Isotropy at fourth order in k is observed for all DnQm lattice stencils. 
This is also true for PK and KU stencils. Other stencils show anisotropic 
discretisation errors. While one may expect Eq. [13] to be anisotropic, due to 
the simplicity in construction, Eq. [15] and [17] use all 27 points of the cubic 
cell, but are still not isotropic at leading order in error. None of these stencils 
provides isotropic error at sixth order. 

In Fig. [21 we plot isocontours of the Laplacian in Fourier space along 
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(d) Shinozaki-Oono, Eq. (c) 'Equally Weighted', Eq. 
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Figure 2: Isocontours of the Laplacian operators in Fourier space, i^(k), at 
two different planes at fc^ = and = ir. Also shown is the isotropic plot 
of k"^ at these two planes for comparison. The color-bar is kept same for 
comparison across the operators. The y-axis is shown only for the first plot 
among each set for clarity. 
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(g) D3Q15, Eq. m 
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(h) Central Difference, Eq. 



Figure 3: Isocontours of the error of Laplacian operators in Fourier space 
defined as -E'(k) = \L(k.) — (— k^)| at two different planes at fc^ = and 
kz = 71 (same as in Fig. |2]). The color-bar is kept same for comparison 
across the operators. The y-axis is shown only for the first among each set 
for clarity. 
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the planes = 0, vr, to visually represent the degrees of anisotropy beyond 
quartic order. Both D3Q27 and Patra-Kartunnen Laplacians exhibit similar 
behavior and are the best of the set. The Laplacian introduced by Kumar also 
exhibits comparable behavior. Both Shinozaki-Oono and 'equally weighted' 
Laplacians are poor approximations at large wavenumbers and as seen from 
Eq. [T^- [TTlthey are not isotropic at quartic order. Both D3Q19 and D3Q15 
are isotropic at quartic order with smaller stencils consisting of 19 and 15 
points. The smallest of central difference stencil is the most anisotropic of 
the set, as seen in Fig. [2hl Laplacians which have identical isotropic errors 
at quartic order may still have different stability properties when employed 
in numerical algorithms. This is part of ongoing work, and will be reported 
elsewhere. Our method naturally generalizes to other derivatives. For in- 
stance, the gradient of ip{r) can be obtained from a Taylor expansion of the 
lattice transform w^iC^^(r + c^). Following exactly the same steps used 
to derive the Laplacian, we see that discretisation error is isotropic. A sys- 
tematic account of gradient, divergence and curl operators obtained by this 



method, as well as the case of genuine anisotropic physics [16|, and higher 



order lattices (17|, ll8[, will be reported elsewhere. 



4. Summary 

Summarizing, we have shown that lattices and weights commonly em- 
ployed in lattice hydrodynamic simulations, provide a computationally effi- 
cient discrete representation of the Laplacian preserving isotropy up to lead- 
ing order error. The weights are derived from the lattice analogue of the 
Maxwell-Boltzmann distribution and are related to Hermite series expansion 
of local Maxwell-Boltzmann distributions. The use of these Laplacians should 
prove beneficial for cell dynamics simulations, hybrid lattice Boltzmann sim- 
ulations, and many other problems, where efficient isotropic discretizations 
of the Laplacian are required. 



Acknowledgement 

One of the authors (SS) wishes to thank the Indian Academy of Sciences 
for financial support and kind hospitality. 



References 

[1] A. Kumar, J. Comput. Phys. 201 (2004) 109-118. 



10 



[2] M. Patra, M. Karttunen, Numer. Meth. Part. Diff. Eq. 22 (2005) 936- 
953. 

[3] Y. H. Qian, D. D'Humires, P. Lallemand, Europhys. Lett 17 (1992) 
479-484. 

[4] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Be- 
yond, Oxford University Press, 2001. 

[5] J. C. Maxwell, V. illustrations of the dynamical theory of gases. part i. 
on the motions and collisions of perfectly elastic spheres, Philo. Mag. 
Series 4 19 (1860) 19-32. 

[6] L. Boltzmann, Wien. Ber. 66 (1871) 275-370. 

[7] S. Ansumah, I. V. Karlin, H. C. Ottinger, Europhys. Lett. 63 (2003) 
798-804. 

[8] 1. Karlin, P. Asinari, Physica A 8 (2010) 1530-1548. 

[9] H. Chen, I. Goldhirsch, S. A. Orzag, J. Sci. Comput. 34 (2008) 87-112. 

[10] M. M. Gupta, J. Kouatchou, Numer. Meth. Part. Diff. Eq. 14 (1998) 
593-606. 

[11] U. Ananthakrishnaiah, R. P. Manohar, J. W. Stephenson, Numer. Meth. 
Part. Diff. Eq. 3 (1987) 229-240. 

[12] W. F. Spotz, G. F. Carey, Numer. Meth. Part. Diff. Eq. 12 (1996) 
235-243. 

[13] A. Shinozaki, Y. Oono, Phys. Rev. E. 48 (1993) 2622-2654. 

[14] V. M. Kendon, M. E. Gates, I. Pagonabarraga, J. G. Desplat, P. Bladon, 
J. Fluid. Mech. 440 (2001) 147-203. 

[15] J. G. Desplat, I. Pagonabarraga, P. Bladon, Gomput. Phys. Gommun. 
134 (2001) 273-290. 

[16] 1. Rasin, W. Miller, S. Succi, J. Gomp. Phys. 206 (2005) 453-462. 

[17] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, F. Toschi, 
Phys. Rev. E 75 (2007) 026702. 

11 



[18] X. Shan, Phys. Rev. E 77 (2008) 066702. 



12 



